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This paper presents an approach to the development of a scaled wind tunnel model for static aeroelastic 
similarity with a full-scale wing model. The full-scale aircraft model is based on the NASA Generic Transport 
Model (GTM) with flexible wing structures referred to as the Elastically Shaped Aircraft Concept (ESAC). 
The baseline stiffness of the ESAC wing represents a conventionally stiff wing model. Static aeroelastic scaling 
is conducted on the stiff wing configuration to develop the wind tunnel model, but additional tailoring is also 
conducted such that the wind tunnel model achieves a 10% wing tip deflection at the wind tunnel test condition. 
An aeroelastic scaling procedure and analysis is conducted, and a sub-scale flexible wind tunnel model based on 
the full-scale’s undeformed jig-shape is developed. Optimization of the flexible wind tunnel model’s undeflected 
twist along the span, or pre-twist or wash-out, is then conducted for the design test condition. The resulting 
wind tunnel model is an aeroelastic model designed for the wind tunnel test condition. 


I. Introduction 

Due to recent strides in the development of light-weight materials, the aircraft industry has been investigating the 
possibility of reducing airframe weight to increase energy efficiency. Reduction of the aircraft weight translates into 
a lower lift requirement and in turn, reduces induced drag, thrust requirements, fuel bum, and cost. These modern 
materials, such as advanced composites, are able to maintain the same load-carrying capacity of conventional air- 
frame material selections. The provided structural rigidity of these materials, however, can be reduced. It becomes 
increasingly important for these modern designs to take into account the aeroelastic interactions of flight aerodynamics 
and the flexible aircraft structures within flight. These aeroelastic interactions can potentially degrade aerodynamic 
efficiency, and thus must be accurately modeled and analyzed. 

A NASA conceptual study titled “Elastically Shaped Future Air Vehicle Concept” was conducted in 2010* to 
investigate the potential benefits of several advanced concepts of a transport aircraft. The study showed that there 
exists potential benefits in shaping wing surface aeroelastic deformation actively in flight with control. In designs 
where stractural flexibility is lessened, active wing shaping control can be used to tailor a wing’s aeroelastic shape. 
The results of the study, however, also showed that conventional flap and slat devices are not aerodynamically ideal as 
control surfaces for active wing shaping control.* 

A novel control surface known as the Variable Camber Continuous Trailing Edge Flap (VCCTEF) system was 
proposed as a new control surface candidate.* Under the Fixed Wing project Active Aeroelastic Shape Control (AASC) 
element, NASA and Boeing are currently conducting a joint study to investigate the application of the VCCTEF 
system^’ ^ on a commercial transport class aircraft. The goal of the VCCTEF study is to investigate the applicability 
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of the VCCTEF as a method to optimize the wing’s spanwise twist shape to establish the best lift-to-drag (L/D) ratio 
during any point within the flight envelope. This offers a significant advantage over the majority of conventional 
commercial aircraft wing designs which are twisted for a set cruise condition and cannot be retailored within flight. 

The VCCTEF is implemented on a model of the GTM^ with structural flexibility of the wing considered, herein 
referred to as the Elastically Shaped Aircraft Concept (ESAC). As investigation of the VCCTEF system continues, 
wind tunnel testing as a method to gauge the potential of the new control surface has been proposed. In the current 
efforts, NASA and Boeing have joined together with the University of Washington Aeronautical Laboratory (UWAL) 
to evaluate the VCCTEF in a subsonic wind tunnel test. A procedure for modeling the development of a wind tunnel 
model configuration using software and numerical tools is developed in order to facilitate decision making with regard 
to wind tunnel testing. 

This paper describes the approach for analyzing and scaling the full-scale ESAC wing to model a wind tunnel 
model configuration prior to inclusion of the VCCTEF. A static aeroelastic model is developed based on the ESAC 
wing’s jig-shape for the candidate wind tunnel model. The model is based on mimicking the aeroelastic behavior of 
full-scale ESAC wing, but higher wing tip deflection is desired of approximately 10% of the wing semi-span. The 
model utilizes a one-dimensional structural model of the the wing structure as a beam in coupled bending-torsion.^’® 
Finite-element method (FEM) will be used to formulate a discretized weak-form solution to the structural equations.®"^ 
Vortex-lattice will be used to conduct the aerodynamic modeling and determine the loads over the wing surface. FEM 
and vortex-lattice are coupled together in structural-aerodynamic loops to generate the aeroelastic model. 

Design of the wind tunnel model configuration is completed with optimization of the wing twist of the undeformed 
wind tunnel model, or wing pre-twist, along the span. This optimization is conducted so that the wind tunnel model will 
experience maximum L/D or minimum induced drag at the wind tunnel test condition. A gradient-based optimization 
using the conjugate directions method and one-dimensional line searches is applied. The resulting wind tunnel model 
is thus a scaled static aeroelastic clean wing model designed for the wind tunnel conditions. 

II. Elastically Shaped Aircraft Concept 

The elastically shaped aircraft concept (ESAC) is modeled as a notional single-aisle, mid-size, 200-passenger 
aircraft based on the NASA Generic Transport Model^ configuration. The GTM is a research platform that includes 
a wind tunnel model and a remotely piloted vehicle, and the geometry of the ESAC is obtained by scaling up the 
GTM wind tunnel model geometry by a scale of 200: 11. Figure 1 is an illustration of the GTM planform. The reason 
for selecting the GTM as a starting point is that an extensive wind tunnel aerodynamic database^ exists that can be 
used in analysis. The benchmark configuration represents one of the most common types of transport aircraft in the 
commercial aviation section that provides short-to-medium range passenger carrying capacities. 



Figure 1. GTM Planform 

In the aeroelastic model of the ESAC, the wing is allowed to freely deform based on reference B757 wing stiffness 
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values and the GTM jig-shape planfom. 

As an active wing shaping control surface, the VCCTEF is implemented on the ESAC jig-shape wing. The VC- 
CTEE is divided into 14 sections attached to the outer wing and three sections attached to the inner wing for each 
side, as shown in Eig. 2. Each 24-inch section has three camber flap segments that can be individually commanded, as 
shown in Eig. 3. These camber flaps are joined to the neighboring sections by using a flexible and supported material 
(shown in blue), which deforms and provides smooth transitions between flap sections without drag producing gaps 
present on most control surfaces of existing aircraft. The flexible skin materials that cover the spanwise camber flap 
sections also constrain the flap deflections such that the relative flap deflections between any two adjacent spanwise 
flap sections are limited. More information on the VCCTEF concept can be found in references.^’ ^ 



Figure 2. ESAC with VCCTEF 


Variable Camber Flap 



Figure 3. Variable Camber Flap 


A. UWAL Wind Tunnel Test 

By early 2013, UWAL had begun investigation into the construction of a wind tunnel wing concept equipped with 
the VCCTEF configuration. This tunnel test was aimed to analyze the behavior of a highly flexible sub-scale model 
and gauge the usage of the VCCTEF as a control surface. A subsonic wind tunnel test at the UWAL facilities was 
planned. A notional diagram representing how the proposed test will be conducted is shown in Fig. 4. The proposed 
wind tunnel test is a floor-mounted test where the semi-span wing model is limited based on the height of the wind 
tunnel test chamber at approximately six feet. For a model of this size, the wind tunnel dynamic pressure is chosen to 
be ^„,„ = 20||. 
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Figure 4. UWAL Wind TUnnel Test Concept 

As joint efforts between NASA, Boeing, and UWAL proceeded, candidate wind tunnel planforms needed to be 
designed. A preliminary wind tunnel model concept with five sections of VCCTEF is represented in Fig. 5. Because 
construction of the 17 VCCTFF sections (14 on outer wing, three on inner wing) would be intensive on a sub-scale 
model, a representation of the VCCTFF would be used instead. Instead of construction of 17 segments, the number is 
reduced and the inner high lift sections are not included. 
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Figure S. UWAL Wind Tlinnel Model Concept 
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B. Wing Alone Model 

As a response to the need to design a wing model to meet the UWAL wind tunnel test requirements, an aeroelastic 
model of a full-scale “wing alone” design is initially generated. Because a candidate wind tunnel model configuration 
does not possess any airframe structures other than the wing, a symmetric model based on GTM geometry that has 
no fuselage, tails, engine nacelles, or pylons is created. The symmetry of the model is only maintained to facilitate 
aerodynamic modeling. 

A water-tight geometry is generated by creating a new airfoil section at the fuselage-wing intersection station and 
shifting the jig-shape wing root to the aircraft centerline. The average location of the intersection between the fuselage 
wing box and the wing on the ESAC is estimated to be 6.1708 ft along the Body Butt Line (BBL) from the aircraft 
centerline. The dihedral of the ESAC wing is preserved, and the new wing alone model still possesses a E = 5° 
dihedral. 



Figure 6. Vortex-Lattice Model of the Wing Alone Model 

As an idealization, the curved wing tip of the original jig-shape geometry is also removed and replaced with an 
“idealized straight wing tip”. This is not expected to have much impact on the overall wing alone model’s aerodynam- 
ics and is conducted only to prevent any issues with the vortex-lattice modeling and analysis. 




Figure 7. Wing Alone Model Idealized Straight Wing Tip 

The wing alone reference area is obtained by integrating the chord over the semi-span and multiplying by a factor 
of 2. Let jB represent the coordinate on the wing along the aircraft pitch axis running from the aircraft centerline 
outwards towards the right wing tip (seen in Eig. 8 as the b 2 direction). The subscript / is used to refer to quantities 
related to the full-scale wing alone model. 

bf 

Sref.f = 2 r c{yB)dyB « 1640.8 ft^ (1) 

Jo 

In the wing alone model, the wing reference area does not include the “fictitious” area that would have been 
covered by the fuselage of the full-scale wing generally included in trapezoidal area estimate of a wing. 

The mean aerodynamic chord is also obtained through integration. 

*’/ 

Cf = ^ c\yB)dyB « 17.0991 ft (2) 

Sref.f Jo 

The wing aspect ratio is determined using the span and reference area. 

bj 

ARf= = 7.6000 (3) 

Sref.f 
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The taper ratio is obtained using the root and tip chord values. 



5.5990 ft 
28.7122 ft 


0.1950 


( 4 ) 


III. Wing Structural Modeling 

A structural model of the wing using beam theory is developed which is later incorporated into a fully coupled 
structural-aerodynamic aeroelasticity model. 


A. Reference Frames 


A 



Figure 8. Aircraft Reference Frames 


Figure 8 illustrates three orthogonal views for a typical aircraft and several associated reference frames. These refer- 
ence frames are useful in developing the structural models of the lifting surfaces of an aircraft, although the coordinate 
frames associated with the aircraft wings are primarily used in this analysis. The aircraft body-hxed reference frame B 
is dehned by the unit vectors b\, b 2 , and b^ which are aligned with the aircraft roll, pitch, and yaw axes, respectively. 
The reference frame C is aligned with the right wing’s elastic axis and is dehned by the unit vectors ci, C 2 , and C 3 . 

Let A be the sweep of the elastic axis. The B frame can be related to C through three successive rotations: 1) the hrst 

/ 

rotation about bj, by an angle of | + A to generate an intermediate reference frame B dehned by the unit vectors bi, 

f / / 

b 2 , and b^ (not shown), 2 ) the second rotation about £>2 by the dihedral angle T of the elastic axis that results in the 
intermediate reference frame C' dehned by the unit vectors Cj, C 2 , and C 3 (not shown), and 3) the third rotation about 
Cj by an angle of n to result in the reference frame C. The transformation can be represented by a series of coordinate 
rotations expressed as 



bi 


— sin A —cos A 0 


cosT 0 

sinT 


b2 

= 

cos A —sin A 0 


0 1 

0 


. *3 . 


0 0 1 


-sinT 0 

cosT 



— sin A cos r cos A 

sinAsinT 


Cl 


= 

cos A cos r sin A 

- 

cos A sinT 


C2 



-sinT 0 


— cosT 


. ^3 . 

The analysis can 

be repeated for the left wing. 

The reference frame L 


1 0 

0 -1 

0 0 


0 

0 

-1 



Cl 


C2 


. ^3 . 


( 5 ) 


and is dehned by the unit vectors di, d 2 , and d^. The B frame can be related to D through three successive rotations: 
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1) the first rotation about —bj, by an angle of | + A to generate an intermediate reference frame B” defined by the 

// // /f n 

unit vectors by, b 2 , and b^ (not shown), 2) the second rotation about by the dihedral angle F of the elastic axis that 

III 

results in the intermediate reference frame D defined by the unit vectors d j , d 2 , and d^ (not shown), and 3) the third 

I 

rotation about by an angle of 7t to result in the reference frame D. The relationship can be expressed as 


bi 


— sinA cos A 0 


cosF 0 sinF 

bi 

= 

— cos A —sin A 0 


0 1 0 

. *3 . 


0 0 1 


— sinF 0 cosF 


sinAcosF 

— cos A 

sinAsinF 


di 

cosAcosF 

sinA 

cos A sinF 


di 

-sinF 

0 

— cosF 


d3 


1 0 

0 -1 

0 0 



di 


d2 


d3 


(6) 


B. Elastic Axis 

An analysis of the combined motion of the left wing is conducted in the present section, and the motion of the right 
wing is considered to be equivalent for symmetric flight. This analysis is equivalent to that in a previous study and is 
included for completeness.® 

Let X represent the coordinate along the elastic axis of a wing running from root to tip. The wing pre-twist angle 
7 (x) thus represents the incidence of the airfoil section at the corresponding elastic axis coordinate. A typical wing pre- 
twist varies from nose-up at the wing root to nose-down at the wing tip and is commonly referred to as a “wash-out” 
twist distribution. 

The internal structure of a wing is typically composed of a complex arrangement of load carrying spars and wing 
boxes that carry the stresses and strains introduced by aerodynamic forces and aeroelastic deflections. For this analysis, 
an equivalent beam approach is used which models the wing’s elastic behavior using equivalent stiffness properties. It 
is a common approach in analyzing aeroelastic deflections^ and can be used to analyze high aspect ratio wings with 
good accuracy. The effect of wing curvature is ignored and straight beam theory is used to model the wing deflection. 
The axial or extensional deflection of a wing is also generally very small and is neglected. 



Figure 9. Left Wing Reference Frame 

Consider an airfoil section on the left wing as shown in Fig. 9 undergoing bending and torsional deflections. Let 
(x,y,z) be the coordinates of point Q on the wing airfoil section. Then the undeformed local airfoil coordinates of 
point Q are 

y 

z 

where rj and ^ are the local airfoil coordinates 


cosy —sin 7 
sin 7 cos 7 


1 




(7) 
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Differentiating with respect to x gives 


7 x 

1 

= 7 

— sin 7 

— cos 7 

J? 


-ZY 


cos 7 

— sin 7 



yy 


( 8 ) 


Let © be a torsional twist angle about the x-axis, positive nose-down, and let W and V be flapwise and chordwise 
bending deflections of point Q, respectively. Then, the rotation angle due to the elastic deformation can be expressed 
as 

^{x,t) = &di-Wjtd2 + V^d3 (9) 

where the subscript x denotes the partial derivatives of 0, W, and V. 

Let (xi ,yi ,zi) be the coordinates of point Q on the airfoil in the reference frame D. Then the coordinates (xi,yi,zi) 
are computed using the small angle approximation as 


Xl(x,f) 


X 


4> X {yd2 + zd2,).di 


1 

1 

1 

1 

yi(x,f) 

= 

y + F 

+ 

X {yd2 + zd2,)-d2 

= 

y + V -z@ 

. zi{x,t) 


z + w 


_ (j> X (yd2 + zd3).d3 


z + W -by© 


Differentiating xi, yi, and zi with respect to x yields 


X\,x 


1 - yVxx + zyVx - zWxx - yY'^x 

yi.x 

= 

-ZY + yx-z©x-yY® 



YY +VYx+y©x-zY® 


Neglecting the transverse shear effect, the longitudinal strain is computed as" 

dsi—ds 

ds Sx 

where 

= \/i+yi+2? = \/i + (/+z2)(7)^ 


(11) 


( 12 ) 

(13) 


s\a = ^Jx\^x+yu+A,x 

= \J si - 2yVxx - 2zWxx + 2(y2 + z^)i@x + (-rpj; - 1)^ + (yi,x + zjY + -yjY 

Ignoring the second-order terms and using Taylor series expansion, is approximated as 

, -yVxx-zWxx+{f + Z^)Y®x 
ilx ~ 

The longitudinal strain is then obtained as 

-yVxx - zWxx + (y^+z^)r'©x 


e = • 


-y 


Si 

..2 I 


i + (y^+z^)(Yr 


Vxx-Z 


l + (/ + z')(7)' W„ + (/+z2)yy l + (/+^2)(y) 


©V 


For a small wing twist angle 7 , (7 ~ 0. Then longitudinal strain can be expressed as 

e = -yFxj - zWxx + ( 7 ^ + ) 7 ©;c 

The moments acting on the wing are then obtained as * ^ 


dydz 


■ Mx ' 


GJ@x 


{y^+z^){Y +©x) 

My 

= 

0 

—z 

. Mz _ 


0 

-y 


GJ + EBi{Yf -EB2Y -EB2Y 


-EB2Y 

—EB^Y 


EL 


yy 


-EL 


yz 


-EL 


yz 


Elz: 



' ©X ' 


^xx 


. . 


(14) 


(15) 


(16) 


( 17 ) 
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where E is the Young’s modulus; G is the shear modulus; 7 is the derivative of the wing pre-twist angle; lyy, ly^, and 
Iyy_ are the section area moments of inertia about the flapwise axis; J is the torsional constant; and B\, B2, and B3 are 
the bending-torsion coupling constants which are defined as 



The strain analysis shows that, for a pre-twisted wing, the bending deflections are coupled to the torsional deflection 
via the slope of the wing pre-twist angle. This coupling can be significant if the wash-out slope 7 is dominant as in 
highly twisted wings such as turbomachinery blades. 

C. Aeroelastic Angle of Attack 

The aeroelastic angle of attack is the effective angle of attack of a flexible wing section that is undergoing aeroelastic 
deformation defined by elastic axis twist 0 , flapwise bending W, and chordwise bending V. It can be calculated by 
solving for the relative velocity of air as it approaches a wing section perpendicular to the elastic axis. The aeroelastic 
angle of attack encompasses a wing section’s rigid local angle of attack and the contribution due to wing elastic 
deformation, and also governs the aerodynamic forces and moments on a local wing section. 

The local angle of attack depends on the relative approaching air velocity, the rotation angle (j>, and the relative air 
velocity in turn also depends on the deflection-induced velocity. The velocity at point Q due to the aircraft velocity 
and angular velocity in the reference frame D is computed as 

vq = v+cox r= {ubi+vb2 + wb 3 ) + {pbi+qb2 + rb 3 ) X {-Xabi-yab2-Zah) 

= {u + rya-qza)bi + {v-rxa+pza)b2 + (w + qxa-pya)b3 ( 19 ) 

where (xa,ya,Za) are the coordinates of point Q in the aircraft body B frame relative to the aircraft center of gravity 
(C.G.) such that is positive when point Q is aft of the aircraft C.G., ya is positive when point Q is towards the left 
wing from the aircraft C.G., and Za is positive when point Q is above the aircraft C.G. The aircraft velocity is (m, v, w) 
in the aircraft body axes, and (p,q,r) are the aircraft angular velocity components. 

This can be expressed in the left wing frame D as vq = x,di +ytd2 + Ztdj, the local velocity due to aircraft 
rigid-body dynamics. 


Xt 


— {u + rya — qza) sinAcosF— {v — rxa+ pza) cos AcosF— {w + qxa — py a) sinF 

yt 

= 

~{u + rya - qZa) cos A-|- (v - rx„ + pza) sin A 



{u + rya — qZa) sinAsinF -f {v — rxa + pZa) cos AsinF— {w + qxa~ pya)cosT 


For a trim case where j 5 =Q, p = q = r = Q, then 


X, 


— MsinAcosF— wsinF 

yt 

= 

— Mcos A 

. 


MsinAsinF — wcosF 


( 21 ) 


The local velocity at point Q due to both aircraft rigid-body dynamics and aeroelastic deflections in the reference 
frame D is obtained as 


V = vq + ^ X p = v.,di+ Vyd2 + v.d-i 


d\ d2 di 


x,-{z + W+ y®)Wx, -{y + V - z&)Vxt 
yt + Vt- {yVx + zWx)Vx, - (z + IT +y©) 0 , 
Zt + W,- {yVx + zWx)Wx, + {y + V- z&)&, 


( 22 ) 


where (x,y,z) are the coordinates for the point Q in the reference frame D without any aeroelastic deflection. 

For static aeroelasticity, all the velocity components of the aeroelastic deflections are set to zero. Thus, Vx = Xt, 
Vy = yt, and = Zf 

In order to compute the aeroelastic forces and moments, the velocity must be transformed from the reference 
frame D to the airfoil local coordinate reference frame defined by The transformation can be performed 
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using successive rotation matrix multiplication operations as 


Vt] 


cos Vx sin Vx 0 

— sinVj; cosVt 0 
0 0 1 


coslVjc 0 sinlV^ 
0 1 0 
— sinWv 0 cosWv 


1 0 0 
0 cos ( 0 + 7 ) sin(© + 7 ) 

0 — sin(©+ 7 ) cos(©+ 7 ) 

cos Vx(vx cos "Wx + 17 sin W«) + Vy sin V, 
cos ( © + 7 ) [— sin Vv ( v.v cos Wt + Vj sin Wr) + Vy cos 14 ] + sin ( © + 7 ) ( — Vj: sin IVv + t'z cos Wr ) 

— sin(© + 7 ) [— sin 14 (v^: cos W 4 + sin V 4 .) + Vy cos I 4 ] + cos (© + 7 ) ( — V;^ sin "Wx + cos H 4 )z; 

Vx + Vyl4 + V^x 

— Vx [14 + H4(© + 7 )] + Vy + Vj [(© + 7 ) ~ VxWt] 

Vx [—Wx + Vx{&+Y)] — t'y(©+ 7 ) + Vz [1 + VxWx{&+ Y)] 



X, 


yt 

- 


cos 114 ) 


( 23 ) 


for small deflections. 

The local aeroelastic angle of attack on the airfoil section due to the velocity components Vrj and in the 
reference frame D is computed as 


Ur = 


+ Wi + Av^ + Wi + Wi {v^+ Wi)Avri 


V-n + 


n 


where w, is the downwash due to the three-dimensional lift distribution over a finite-aspect ratio wing. 
The velocity components are 

= MsinAsinT — wcosT 
Av^ = Vx [—Wx -f 14 (© -f 7 )] — Vy (© -f 7 ) 


Vt] = — mcosA 


^V-q — —Vx [14 + H4 (© -f 7 )] + Vz [(© + 7 ) ~ yxWx] 

The local aeroelastic angle of attack is evaluated as 

MsinAsinT — wcosT + (— usinAcosT — wsinT) [— H4-I-14 (©-f 7 )] + mcos A(©-|- 7 ) -fi 


(24) 


(25) 

(26) 

(27) 

(28) 


«£ = -- 


MsinAsinT — wcosT -f w,- 


M^COS 


2A 


u cos A 

— (— MsinAcosT — wsinT) [ 14 - 1-114 (© + 7 )] 


-I- (MsinAsinT— vvcosT) [(©- 77 ) — V4H4 


(29) 


Assuming a trim case, let u « 14, and w « 14,a. Neglecting chordwise bending components, V, 14, and also 
neglecting the three-dimensional finite-wing effect, w,, allows Uc be expressed as 


0!c = - 


sinAsinT acosT sin A cos T -7 a sin T 


cos A cos A 

sinAsinT— acosT 


cos A 


H4-©-7 


(MsinAcosT-7 wsinT) ( 114 ©- 71147 ) -7 (MsinAsinT — wcosT) (©-7 7 ) > (30) 


ac = -- 


MCOs2 A 

sinAsinT acosT sinAcosT -7 asinT 


cos A 


cos A 


cos A 


H4-©-7 


' — sin^ AsinTcosT- 7 a sinAcos^T — asinAsinT2-7a^sinTcosT', , 

+ 1 ) (H4© + H47) 


COs2 A 


COs2 A 


-7 


— sin^ Asin^T -7 asinAsinTcosT a sinAsinT cosT— a^cos^T 


cos^A cos^A 

Eliminating higher order terms results in the aeroelastic angle of attack expressed as 

cosT 

Uc = — 7 — tanAsinT -7a © — H4tanAcosT 

cos A 


(©-77) 


(31) 


( 32 ) 
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This can be re-expressed after applying small angle approximation in terms of partial derivatives as 


ccc{x,y,z) 


dUc 

In 


dUc 

da 




dOc. 

dW. 




dUc 


© 


dUc 

d\ 


—7— tan AT 


dac 

da 

dac 


1 

cos A 


= — tan A 


dac _ 
d@ ~ 

The aeroelastic deflections terms Wx and © contribute to aerodynamic stiffness. 


(33) 

(34) 

(35) 

(36) 

(37) 


D. Coupled Bending-Torsion Equations 


Without considering chordwise bending of the wing, the equilibrium conditions for bending and torsion are expressed 
11 

= -mx (38) 


as 


dx 


d^My 

dx 


- f _ 

— Jz z\ 

ox 


(39) 


where nix is the pitching moment per unit span about the elastic axis, is the lift force per unit span, and niy is the 
bending moment per unit span about the flapwise axis of the wing. 

Because the structural modeling is intended for use in a static aeroelasticity model, a steady-state aerodynamics 
model is used. Aerodynamic information can be obtained through vortex-lattice modeling to develop the forces and 
moments for coupled bending-torsion of a flexible wing. 



Figure 10. Airfoil Forces and Moments 

Neglecting the effect of downwash that is caused due to lift generation over a three-dimensional finite-wing, the 
lift coefficient over the span of a clean wing assuming linear aerodynamics is as follows; 

Cl{x)= CLa (.r ) ( a ) (40) 

where is the aeroelastic angle of attack as shown in Fig. 10, assumed to be constant for airfoil cross sections 
perpendicular to the elastic axis and only a function along the wing such that = ac{x). 

The aeroelastic angle of attack can be expressed as contributions due to aircraft rigid-body angle of attack a as 
well as the contribution due to aeroelastic deformation. Let ar represent the contribution the aeroelastic angle of attack 
due to rigid-body considerations including airfoil shape, and represent the effect on the local aeroelastic angle of 
attack due to aeroelastic deformation at the aerodynamic center of the airfoil section. Based on Eq. 33 and neglecting 
chordwise bending and assuming dihedral is small, ar and can be represented as 

ar{x) = 
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ac{x) = 


ar{x) + ae{x) 
cos A 


(42) 

(43) 


where both a^. and are about the pitch axis direction, positive nose-up. 

The rigid and elastic lift coefficient contributions to the local sectional lift coefficient of the elastic axis airfoil cross 
sections can be expressed as 

cl,(x) =cl„(x)^^ (44) 


cl,(x) =cl„(x) 


cos A 
Oie{x) 
COS A 


(45) 


It is also important to note that the elastic contribution to the local aeroelastic angle of attack can be repre- 
sented based on the partial derivatives calculated in Eqs. 36 and 37. Given a deformation characterized by elastic axis 
twist © and vertical bending slope Wx, the elastic contribution to the aeroelastic angle of attack can be calculated as 


ae{x) = — 0(x)cosA — Wt(x) sin A 


(46) 


where Ue is about the aircraft pitch axis. This applies for the static case using the assumptions and simplifications 
applied in derivation of Uc- 

The steady-state drag coefficient can be modeled by a parabolic drag polar as 


cd{x) = cdo(x) -|-k(x)ci(x) 

where is the section parasitic drag coefficient and k is the section drag polar parameter. 
Likewise, the pitching moment coefficient about the aircraft pitch axis can be represented as 

c„,(x) = C„,„^(x) -b ^Cz,(x) cos A 
c(x) 


(47) 


(48) 


where e is the location of the aerodynamic center relative to the elastic axis along the body axis, positive when the 
aerodynamic center is forward of the elastic axis, and is defined about the pitch axis, positive nose-up. 

Expanding Eq. 48 using the aeroelastic angle of attack Ue definition in Eq. 33 produces 


/X /X , .fdUc dUc dac,,\ 

= C..J.) + ®+ 

which allows us to define the following quantity 

, , e(x) , , / dUc dUc \ 

The lift force, drag force, and pitching moment about the aircraft pitch axis are expressed as 

/ = ciq^o cos Ac 

d = cdQ'ooCOS Ac 
m = c,nqo=c^ 


(49) 

(50) 

(51) 

(52) 

(53) 


where cos A takes into account the correction due to the elastic axis sweep but is not needed in the pitch moment 
calculation since Eq. 48 is already about the pitch axis. 

The forces and moments in the local coordinate reference frame are obtained as 


/" = (/cos a -|-c/sina)r + {dcosa — I sin a) sin A 

/" = (c/ cos a — /sin a) cos A 
/" = / cos a -be/ sin a — (c/cosa — / sin a) sin AT 
ml = — mcosA 


(54) 

(55) 

(56) 

(57) 
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m" = msinA 
= mcos Ar 


(58) 

(59) 


For a model with only flapwise bending and torsion considered, the beam deflection analysis the only aerodynamic 
force and moment terms that have an effect are the terms /", m“, and m". For this analysis, the aerodynamic force and 
moment terms are thus considered to be 


/" « ctqooCOS^Ac 

m“ Ri — Cm^ooCOS^ Ac^ 

dc„, 2 

— (^oo sin A cos Ac 
ox ox 


(60) 

( 61 ) 

(62) 


where an additional cos A term considers the change in the direction of q^o over the wing due to sweep. 

Inserting Eq. 17 and the force and moment terms Eqs. 60-62 into the governing equilibrium equations Eqs. 38 and 
39, the following equations can be used to describe the coupled bending and torsion motion of the wing: 


dx^ 


{-EB2Y&x + EIyyW„) = -mW, 


2 dCj^j 2 2 

- mecg&tt + Ciqoo cos Ac — q^o tan Acos Ac 

* rlr 


(63) 


l{[GJ + EBarf 


&X - EB2YWxx 1 = mrl&tt - mecgW,, + c,nq< 


,cos^Ac^ 


(64) 


Although the wing alone model is modeled as a symmetric wing in a horizontal plane, the actual wing tunnel model 
will be only a semi-span wing mounted in a vertical plane. Gravitational forces on the wing alone and candidate wind 
tunnel model are thus ignored. There are also no engines on the model, and thus the only forces and moments being 
considered in the aeroelastic model are from aerodynamic and inertial sources. 


IV. Finite-Element Modeling 


The development of the coupled bending-torsion partial differential equations describing the wing allows for wing 
bending and torsional deflections to be solved. EEM^ is used as a numerical technique that uses locally-defined basis 
functions to numerically approximate the solution of the governing partial differential equations. The EEM is used to 
discretize the wing structure into n equally spaced one-dimensional elements. The bending and torsional deflections 
can be approximated as 

@{x,t) = '£@i{x,t) (65) 

/=1 
n 

W{x,t) = Y^Wi{x,t) (66) 

(=1 

where i refers to the /-th element. 

Eor each element, the bending and torsional deflections are approximated as 


©,(x,f) = I/, '01,(0 + V2{x}92i(t) 


V/l(x) \j/2(x) 


01,(0 

02,(0 


= Ng{x)9i{t) 


(67) 


Wi{x,t) = 


(j)l{x)wi.{t) + ^ 2 (x)w.(t) + (j) 2 {x)w 2 ,{t) + 04(x)w2,(f) 


^l(x) 02 (-v) h{x) 04 (jc) 


1V1,(0 
Vl’l,(0 
W2,(f) 
W^2,(0 . 


= N»,{x)wi{t) 


( 68 ) 


where the subscripts 1 and 2 denote values at nodes 1 and 2, and i//(x), j = 1,2 and 0^(x), k = 1,2, 3, 4 are the linear 
and Hermite polynomial shape functions 

y/j (x) = 1 - y (69) 
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(70) 


V/2(jc) = 


0i(x) = 1-3(^)2 + 2(^)3 


02(x) = /^y-2(^)2 + (^)3 


03(x) = 3(y)2-2(-)3 


(71) 

(72) 

(73) 

(74) 


where x S [0, 1] is the local coordinate and I = ^ is the element length. 

The weak-form integral expressions of the coupled bending-torsion partial differential equations are obtained by 
multiplying the equations by Ng (x) and A(J(x) and then integrating over the wing span. The aerodynamic coefficients 
are expanded here based on the aeroelastic angle of attack representation in Eq. 33 and using Eq. 49. This yields 


f-Jo ^ + N'eet-EB27Nlwt}dx 

n pi n pi 

^ / NlmrlNQdi-mecgN„Wi)dx+Y^ / N^ 

i=\Jo i=\Jo 


c„,o -f - cos Aci„ I — iVe 0,- -f 


qoocos^ Kc^dx (75) 


E / Nl'^{-EB 27 Ngei + EIy,Nlwi)dx = 
i=QjQ “ 

n n pi ' 

Y^Nl{pAN„Wi + pAecgNeBi)dx+Y^ / cl„ 

i=\ i=i 7 o 

-u 


dUc , <3«c , „ , dUc' 


qooCOS^ Acdx 


N — 

x ’ w; 7 

ax 


c,„o -f - cos Acl„ I -^NgOi + ^N„Wi 


tanA^oo cos^Ac^dx 


(76) 


The expressions of the left hand sides can be integrated by parts upon enforcing the boundary conditions, resulting 


/' 


lvT_(_£B27A^e0(+£7v>7Vw 


Z’/ 

i7,)£/x= / 

Jo 


-EB27NQ 6i + EIyyN„Wi)dx 


The elemental mass matrix, stiffness matrices, and force vector are then established as 


Mv- = / m 


/' 


rlNlNe -CcgN^N„ 
-ecgNlNg NlN^ 


dx 


K.= 


f 


GJ + EBi{7f 


1 f'T fl ~ 

-EB 27 N^N„ 

Jo 


/ //t- 

-EB27NJ 

K 

El n''^ n” 


Ka,= 


ecL„ cos A^NjNe ccl„ cos A 


. ^ aA a / 


— c, 

^La ^0 




qooCOS^ Acdx 




Er . = / Q'ooCOS Ac 


^MqCN Q 

ClK 


0 

-CmoCtanAA^7 


(Jx 


-f g'ooCos'^ Ac 


-CmgCtanAA^ 


(78) 


(79) 


(80) 


(81) 


(82) 
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where is the structural stiffness matrix and Ka is the stiffness matrix due to the result of aerodynamics. The moment 
component due to lift eccentricity that contributes to bending terms in Ka are neglected. 

The globally assembled system is described by the matrix equation 


T K^Xg — Fy KaXg 


(83) 


where Xe = 


tT 


01 W1 Wi 02 W2 W7 


0, 


Wi W: 


0«+i w„+i w, 


n+1 


Equation 83 represents the governing equation for solving the structural deflection of a flexible wing given aero- 
dynamic force and moment inputs. By setting Xe = 0, the equilibrium solution forxg can be obtained through inverting 
the stiffness matrix and pre-multiplying the force matrix. Information can then be extracted from the solution including 
the wing’s deflection along the elastic axis of the wing © and W. It can also be used to calculate the elastic contribution 
to the aeroelastic angle of attack in Eq. 33. 

Without considering the effect of aeroelastic coupling that results in the aerodynamic stiffness matrix Ka, a wing 
structural deflection can be solved by 


Xe=K-^Fy 


(84) 


Aeroelastic deflection of the wing can be calculated by including the aerodynamic stiffness. The aeroelastic solu- 
tion is represented by the term Xe to differentiate it from Xg, which is structural deflection calculated without using the 
aerodynamic stiffness matrix. 


Xe={Ks + Ka)-^Fr 


(85) 


The term Fy represents the force matrix in the EEM that is constructed solely based on the aerodynamic character- 
istics of the undeformed wing. Thus, Eq. 85 represents a method in which the aeroelastic deformation of a wing can 
be solved solely using the rigid wing aerodynamic loads and properties. While this is a powerful framework, another 
approach exists that involves updating the force matrix in the EEM instead of utilizing the rigid wing properties. 

Let F represent the load on the wing due to the rigid planform and the incremental effect due to aeroelastic 
deformation. The force F can be represented as 


F = Fy-KaXe ( 86 ) 

where Xg is the deflection of the wing. Aeroelastic deformations can be calculated by an iterative technique where 

= K^^F^, Xg^* — > Xg as k increases (87) 

where the force vector F^ is estimated using computational aerodynamic modeling tools to calculate the loads on 
the deformed geometry. This approach, which is used in this study, requires an aerodynamic modeling tool to be 
run coupled with the EEM model. One of the advantages of using this method is that the aircraft and sectional 
characteristics are estimated for the deformed geometry, instead of using the assumption that they remain constant at 
the rigid wing values. This can also provide a better estimate of the effect of Ka than the analytical model developed. 

V. Vortex-Lattice Aerodynamic Modeling 

Vorview is a computational tool used for aerodynamic modeling of aircraft configurations using vortex-lattice 
method.*^ Based on lifting line aerodynamic theory, Vorview provides a rapid method for estimating aerodynamic 
force and moment coefficients. Geometric input vehicle configuration are constructed within Vorview by discretizing 
the surface into a series of panels, which are then represented by placement of spanwise and chordwise locations 
of bound or horseshoe vortices. Vorview computes the vehicle aerodynamics in both the longitudinal and lateral 
directions independently, and these can be combined to produce the overall aerodynamic characteristics of the vehicle 
at any arbitrary angle of attack and angle of sideslip. 

Vorview is considered a medium fidelity tool, and limitations associated with vortex-lattice modeling in general 
apply to Vorview aerodynamic analysis. The drag prediction by Vorview is most reliable only for induced drag 
prediction due to the inviscid nature of any vortex-lattice method. Prediction of viscous drag due to boundary layer 
separation and wave drag due to shock-induced boundary layer separation are generally not conducted by vortex- 
lattice, and viscous drag must be estimated using other methods. 

In addition to force and moment analysis, Vorview can provide a rapid estimation of aerodynamic derivatives 
including dynamic derivatives due to angular rates. These aerodynamic stability and control derivatives are useful 
in analyzing the stability and handling characteristics of an aircraft configuration. Owing to the computationally 
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efficient vortex-lattice method, aerodynamic derivatives can be estimated in Vorview fairly quickly. A flight dynamic 
model for a given vehicle configuration can be easily developed with Vorview using the results from these stability 
and handling analyses. Vorview has been validated by both wind tunnel data"* as well as the NASA Cart3D tool,^^ 
which is a high-fidelity inviscid (Euler) CFD analysis code targeted at analyzing aircraft performance in conceptual 
and preliminary aerodynamic design. In general, both Vorview and Cart3D seem to have similar predictive capabilities 
when compressibility is not a factor. 

Figure 1 1 illustrates an aerodynamic model of the GTM in Vorview. 



Figure 11. GTM Aircraft Model in Vorview 

In this study, Vorview will be utilized as the primary tool for conducting aerodynamic modeling for the aircraft 
configurations. Total aircraft characteristics as well as sectional data along the aircraft wing surfaces can be post- 
processed from Vorview. 


VI. Automated Geometry Modeling Tool 

An automated geometry generation tool is developed in Matlab that is used to close the loop between the structural 
and aerodynamic modeling needed to generate an aeroelastic model. The geometry generation tool uses structural 
deflection data that is computed by the FEM model and applies it to the undeformed aircraft wing geometry to reflect 
static aeroelastic deflections. The vehicle geometry modeler directly outputs a geometry input file that can be read by 
Vorview when computing an aeroelastic solution. 


A j 1 








Figure 12. GTM Coordinate Systems 

Consider the reference frames in Fig. 12. The coordinate reference frame {xA^yA,ZA) defines the Body Station 
(BS), the Body Butt Fine (BBF), and the Body Water Fine (BWF) of the aircraft, respectively. The coordinate 
reference frame {xv^yv:Zv) is the translated coordinate system attached to the nose of the aircraft such that xy = 
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xb — 13.25 ft, yv = yn, and zv = zb~ 15.8333 ft. This reference frame is used by vortex-lattice aerodynamic modeling 
tool. The aircraft body reference frame (xB,yB,ZB) is the same B coordinate system dehned earlier in Fig. 8 by the unit 
vectors bi, £>2, and £>3. The B coordinate frame is attached to the aircraft center of gravity (CG) such that xb = xv—xy, 
yB=yv — yv, and Zb = Zy—Zy, where {xy.yyZy) is the coordinate of the CG in the {xy.yyZy) reference frame. 

The vehicle geometry modeler has access to the outer mold line of the aircraft geometry. It is capable of applying 
geometric transformations onto the outer mold coordinates of the wing’s jig-shape to simulate aeroelastic deflection. 
Neglecting chordwise bending deflection and utilizing the coordinate system of the left wing developed earlier (coor- 
dinate frame D), the aeroelastic deflections in bending and torsion are in expressed in a vector form as 

W_,d2 ( 88 ) 

Ar = —W sinlT^flf 1 -f W cos Wxdi (89) 

The coordinate reference frame {x,y,z) of the left wing is related to the coordinate reference frame {xy,yy,zy) by 
the following relationship 


di ' 


— sin A COS F 

— cos A sin F 

-sinF 


‘ £>i ' 


di 

= 

— cos A 

sin A 

0 


£>2 


d3 


sinAsiny 

cos A sin F 

— cosF 


. *3 . 




— sin A cos F 

— cos A cos F 

-sinF 


-Vi 



= 

— cos A 

sinA 

0 


V2 



sin A sin F 

cos A sin F 

— cosF 


-V3 


where (vi, V2, V3) are the unit vectors for the Vorview coordinate reference frame {xy^yy^zy) ■ 

Thus, the aeroelastic deflections result in a wing twist expressed as an incremental angle of attack Aa (positive 
nose-up), a horizontal deflection Ayy (positive deflection towards wing tip), and a vertical deflection Azy (positive 
displacement upward) as follows: 

Aa = — ©cosAcosF — WtsinA (91) 

Ayy = — VTsinWtCos AcosF — VTcoslTjcCosAsinr (92) 

Az„ = — ITsinlTvSinr-l-lTcoslTj^cosr (93) 

A coordinate transformation to account for wing aeroelastic deflections is performed by rotating a wing section 
about its elastic axis by the incremental angle of attack Aa and then translating the resultant coordinates by the 
horizontal deflection Ayy and the vertical deflection Azy. 

Note that the transformation for Aa is equivalent to the value of a^, the local change in the angle of attack for a 
wing section due to aeroelastic deformation represented by Eq. 46, when dihedral F is small. 

VII. Static Aeroelastic Model 

In a standard static aeroelastic model, it is understood that the modeling effort needs to take into account that 
structural deformations during flight will alter the aircraft aerodynamics, and changing the aerodynamics will thus 
change the structural deformations. This realizes an aeroelastic model where coupling exists between the structural 
modeling and aerodynamic modeling approaches. Previous studies have analytically constructed fully coupled aeroe- 
lastic finite-element models that utilize rigid wing lift-curve slopes as an aerodynamic model. This study employs a 
static aeroelastic model that is constructed by utilizing a structural FEM model coupled with a vortex-lattice solution. 

A static aeroelastic code is developed by utilizing the automated geometry generation modeling tool to close the 
loop between the EEM model and the vortex-lattice model using Eq. 87. Eor a model considering only flapwise 
bending and axial torsion, the aeroelastic deflection can be summarized by the quantities of ©(2c), W{x), Wx{x), the 
aeroelastic elastic axis twist, aeroelastic vertical (flapwise) bending, and aeroelastic vertical bending slope respectively. 
These quantities are emphasized to be aeroelastic deflections, while the terms ©(x), W{x), Wx{x) are considered 
structural deflection terms which may or may not be the aeroelastic solution for a given flight condition. Closing the 
static aeroelastic loop causes the the structural deflections ©(x), W{x), Wx{x) to converge to the aeroelastic solution 
©(x), W{x), Wx(x) as iterations are conducted. The structural and aeroelastic deformations can also be represented by 
the elastic contribution to the aeroelastic angle of attack in Eq. 46, or ae(x) and ag(x). 
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Figure 13. Static Aeroelastic Model Concept 

The static aeroelastic model maps an input desired Cl, Mach number M, and altitude h into the respective static 
aeroelastic deflection (©, W, Wx) and the angle of attack a that the flexible wing aircraft would experience when 
leveled at the desired Cl- The following procedure is followed: 

1. Vortex-lattice modeling is conducted on an input geometry at an input flight condition a, M to determine the 
aircraft total aerodynamic quantities, as well as sectional aerodynamic distributions of cl{x), c,„^{x), k{x), 
CLa{x), and x^dx) or the location of the section aerodynamic centers. 

2. The structural FEM model uses the sectional aerodynamic inputs to calculate the wing’s structural deflection 
0(x) and W{x). 

3. The geometry generation tool converts 0(x) and W (x) into the series of deformations ae{x), Ayy (x), and Azv(a), 
and generates a new aircraft geometry with the deformed wing. 

4. A lift curve is generated based on the deflected wing aircraft geometry. The angle of attack a for the value Cl 
is determined and selected for the next iteration. 

5. Steps 1-4 are repeated until Aa between iterations is within a criteria. 

This converged solution is represented by the angle of attack a that the flexible wing model would need to have a lift 
coefficient of Cl- The wing shape at the converged flight condition is the converged aeroelastic deflection ©, W, Wx, 
and cCe- 

The model is used to determine the full-scale wing alone model’s static aeroelastic deflection using the baseline 
stiffness values. A cruise flight condition for the full-scale ESAC is considered to be at Mach = 0.797, altitude 
h = 36,000 ft, with a wing loading of ^ corresponding to a design Cl = 0.510. 

The code is first run restricting any coupled structural-aerodynamic loops. This solution thus represents the case 
where the structural deflections do not affect aerodynamics experienced on the flexible wing, or a model where the 
aerodynamics correspond to the rigid planform only. The deflection results are presented as Wup and @tip, where Wtip 
is the wing tip vertical deflection (positive upwards), and ©„p is the wing tip twist about the elastic axis (positive 
nose-down). These results are summarized in Table 1. 
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Cl 

0.510 

a, deg 

2.279 

W,ip, ft 

2.986 

loo(^), %r 

5.348 

®tip, deg 

-0.351 

/ 2(=>f,p s deg 
1 ft b ft 

-6.286 X 10-3 


Table 1. Structural Deflection Results for Full-Scale Wing Alone Model 

The model is then used to determine the aeroelastic deflection allowing the structural-aerodynamic loops. The 
results are summarized in Table 2. 


Cl 

0.510 

a, deg 

3.064 

%p, ft 

2.740 

lOO(^), % 

4.907 

©rip, deg 

-0.231 

( \ deg 
t ft b ft 

-4.135 X 10-3 


Table 2. Aeroelastic Deflection Results for Full-Scale Wing Alone Model 


VIII. Wind Tunnel Model Scaling 


A candidate wind tunnel model is generated from conducting scaling of the full-scale wing alone model. Scaling 
must be conducted considering several factors and desired characteristics. 


A. Geometric Scaling 

Geometric scaling of the wing alone model is conducted so that the wind tunnel model can fit within the wind tunnel. 
Given a desired wind tunnel model, a geometric scaling factor iiscaie can be determined based on the span of the full- 
scale wing alone model bf and the desired span of the wind tunnel model bg. The subscript w hereinafter refers to the 
sub-scale wind tunnel model characteristics. 

Suppose a wind tunnel height is given to be 6 ft and a desired wind tunnel model semi-span is 5.4219 ft. A 
geometric scaling factor can be determined by: 


bf 2(56.1625 ft) 

flscale 2(5.4219 ft) 


10.3585 


(94) 


The geometric scaling factor can be used to scale all the coordinates of the full-scale wing alone model that is 
used by the geometry generation tool. The reference values for the sub-scale wind tunnel model can also be obtained 
through the scaling factor. 

Sref,w = -^SrefJ = 15.2919 ft^ (95) 

^^scale 

c„,= -^ = 1.6507 ft (96) 

^ scale 

The aspect ratio and the taper ratio remain unchanged. 


AT?* = 7.6000 (97) 

A,, = 0.1950 (98) 

The geometric scaling is not expected to affect the aerodynamics of the models given that all the reference values 
are correctly computed. To verify this, the lift and drag curves for the full-scale model are generated and compared to 
that of the scaled down wind tunnel model. Note that the drag polar only includes vortex-lattice computed drag. 
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Figure 14. Lift and Drag Curve Verification for Full-Scale Wing Alone and Sub-Scale Wind Ttinnel Model, Undeformed and Rigid 


The lift curve and drag polar for the geometrically scaled down wind tunnel model rests almost exactly on top 
of that of the full-scale wing alone model. This indicates that the scaling was done properly such that the total 
aerodynamics of the two models are preserved with geometric scaling. 

B. Aeroelastic Scaling 

Aeroelastic scaling is conducted to determine the scaling factors on the wing’s stiffness properties. From the results in 
Tables 1 and 2, it can be seen that for a flexible wing model, structural-aerodynamic coupling changes the deformation 
result. It can even be modeled that aerodynamic considerations make the wing model stiffer in bending and softer in 
torsion while coupling the two motions,® and this can be deduced from observing the terms in Eq. 81. Because the 
magnitudes of the deformations are generally similar when the total load over the wing is maintained, an uncoupled 
model will be used for preliminary sizing. This approach simplifies the scaling process and allows the torsional and 
bending stiffness of the wing to be analyzed separately due to the fact that the B\ and Bi terms in Eq. 18 are considered 
to be negligible, removing any coupling between the two deformations. 


1. Torsional Stiffness Scaling 

It is known that the structural deformation is calculated using the system equation given by Eq. 83. In the spirit 
of the finite-element analysis previously presented, if a single element modeled with constant structural properties is 
considered with a fixed end, the static deflection equation can approximated by 


e = {ks .torsion) fa.ti 


(99) 


where the value of ks^torsion is related to to the torsional constant J, k^jorsion ^ "'h®re / is the length of the beam 
element. 

In designing the wind tunnel model from the full-scale wing alone, it is desired that the relative elastic axis twist 
of the wing alone and the wind tunnel model are equal. That is, for each beam element: 


If Iw 


Let kj — k^.torsion and ff — fa^torsion^ 


0V. 

If Iw 

jktf) ^ftf ^ {ktff 
If Iw 

If l\v 


(100) 


( 101 ) 

( 102 ) 

( 103 ) 
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(104) 


f‘f _ ftw 
GJf GJ„ 

The force terms are expanded into the aerodynamic contributions, and Se is the reference area for an element based 
on Eq. 82. 

fa, torsion = + Cl - COS A)q'„o COS^ (105) 

The values of the aerodynamic coefficients for the wing alone and wind tunnel model as shown in Fig. 14 are 
equivalent. = Cm^cw — ^niac ^ttd cl,/ = cl.w ~ cl. 


- {cmacCf + Cief cos a) q^COS^ASe,f _ - {Cm^,C„ + CLe»,COsA) qooCOs'^ ASe,» 


GJ 


f 


GJ„ 


( 106 ) 


It is also known that the geometric parameters of the flying wing and the wind tunnel model are related through 
the value Ujca/e thus allowing us to formulate a relationship between the torsional constants. 


- (cm,,, C/ + CLC/ cos A) ^ooCOS^ASgJ _ - (c,„„, [cf/tlscale] +CL [c/AWe] COS A) q^oCOS^ A (Se,f 
“ (cm^^c/ + CLC/ cos A) ^^COS^A^eJ - (c,„„^C/ + ClC/COS A) COS^ A5g J 


GJ 


f 


^•^wnlcale 


GJ„ = 


1 


GJf 


n 


(107) 

(108) 
(109) 


scale 

^ torsion 


Thus the value rtorsion is defined as a scaling factor on the GJ of the flying wing model when scaling for the wind 
tunnel model where rtorsion = ■ 

^ scale 

The static structural deflection for the full-scale and the scaled down model with the application of the scaling 
factor rtorsion are tabulated in Table 3. It is shown that the application of the torsional scaling factor is able to match 
the relative elastic axis twist relative to the span for both models when examined using an uncoupled stmctural- 
aerodynamic model. 


Cl = 0.510 

Full Scale Wing Alone 

Wind Tunnel Model 

a, deg 

2.279 

2.249 

@tip, deg 

- 0.351 

- 0.308 X 10-1 

/ ^^tip \ deg 
\~fT~ ~fT 

- 6.286 X 10-3 

- 5.707 X 10-3 


Table 3. Structural Deflection Results for Wind Tunnel Model with Torsional Stiffness Scaling 

In this analysis, the relative twist per semi-span was preserved in Eq. 100, and this is one of two ways in which 
torsional stiffness scaling can be conducted. Another approach for conducting the scaling would be instead to preserve 
the magnitude of the twist, or 

0/ = 0w (110) 

Following the same derivation process shown before, Eq. 104 becomes 

A/ A ft„iw 


GJf GJ„ 

Substituting in Eq. 105, then 

“ {cntacC f + Chef cos A) qocCOS^ ASejlf _ - {CmacCw + CLewCOsA)qooCOS^ ASe,wk 


( 111 ) 


(112) 


GJf GJ„ 

{^CftiacCf CL6f cos Asj qco cos ASgflf (j'lnac f / ^i-scale\ COsA^^ooCOS ^{^^e,f/^ttcale} h^scale) 


GJ 


f 


GJss 


(113) 
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( 114 ) 


“ {cmacCf + Cte f cos a) qooCOS^ ASejlf - {cm^^c f + CtCf cos A) qooCOS^ ASejlf 


GJf 

If this scaling method is used, the result is 




GJ„ = 


‘‘scale 


GJf 


^ torsion 


(115) 


Thus two possible scaling factors for torsional deflection are developed. For scaling that preserves the relative twist 
relative to the span of the model, a scaling factor rtorsion = can be used. For scaling that preserves the magnitude 

^scale 

of twist, a scaling factor rtorsion = can be used. 

^^scale 


2. Vertical Bending Stijfness Scaling 

Scaling the vertical bending stiffness Elyy is conducted similar to the torsional stiffness scaling. In the development 
for this wind tunnel model, however, the bending stiffness is desired to be scaled such that a 10% tip deflection relative 
to half of the span is achieved. For a beam element undergoing bending, the deflection can be given by 


tv — {l^s, bending) fa, bending 


(116) 




The value of ks bending is given by the relationship kg, bending —jT~ where / is the length of the beam element. 

Let kb = ks^bending ^nd fb = fa, bending- Let an additional scaling factor iibending be defined to increase the tip deflec- 
tion to 10% relative to the model’s semi-span. 


bi-bending , — , 

If Iw 


bibending 


{kbf)-^fbf _ {kbj-'fb.. 


^ bending'^ 


If Iw 

[{EIyy)f/{lf)^] -^fby [{Elyy)„/{1„)^] 


If /h 

fbf{lf)^ _ fb^,{lw)^ 

^bending 


{EIyy)f (EIyy)t, 


(117) 

(118) 

(119) 

( 120 ) 


Vertical bending from aerodynamic sources is due mainly from the sectional lift coefficient. Ignoring the contri- 
bution to the vertical bending force due to the component of pitching moment, the bending force on a beam element 
can be represented as 

( 121 ) 


Iibending 

^bending 


fa, bending — C^^ooCOS A.Se 

CLqooCOS^ ASe,f{lf)^ _ CLqooCOS^ ASe,w{lw)^ 

{EIyy)f [EIyy)„ 

ciqoocos^ ASe,f{lf)^ ciqoocos^ AS e,f {If 


{EIyy)f 
{EIyy)„ = • 


1 


bibendingbi scale 

^bending 


{EIyy)wbi'lcale 

-{EIyy)f 


(122) 

(123) 

(124) 


Thus the value rbending is defined as a scaling factor on the baseline Efy when scaling for the wind tunnel model 
where rbending — 


'‘bending" scale 

Table 4 shows the results for the wind tunnel model’s static structural deflection when computed without any 
structural-aerodynamic coupling and with the application of the scaling factors derived. It can be seen that when the 
scaling factors are applied, a 10% wing tip deflection is achieved on the wind tunnel model. 
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Ci = 0.510 

Lull Scale Wing Alone 

Wind Tunnel Model 

i^^bending ~ 1) 

Wind Tunnel Model 

ij^bending ~ 0.537) 

a, deg 

2.279 

2.249 

2.249 

W,ip, ft 

2.986 

0.289 

0.539 

loo(^), <%r 

5.348 

5.370 

10.000 


Table 4. Structural Deflection Results for Wind Tunnel Model with Vertical Bending Stiffness Scaling 


3. Dynamic Pressure Effects 

A limitation in wind tunnel tests is that wind tunnel facilities may not be equipped to run at Mach numbers as high as 
the cruise condition of a full-scale aircraft. While use of a transonic wind tunnel can be used to test a model at high 
Mach number and high dynamic pressure qoo, operational and usage times at these facilities are very costly. For the 
VCCTEF study, it is not necessary that the wind tunnel model will need to run at high dynamic pressure provided that 
the elastic stiffness is scaled such that the lower dynamic pressure at the wind tunnel test speed is accounted for. 

Let the value q^o.f represent the dynamic pressure of the wing alone model at a cruise Mach number and altitude. 
Let the value represent the dynamic pressure of the wind tunnel model, which is restricted based on the wind 
tunnel test configuration. 

If the wind tunnel dynamic pressure is the same as the wing alone model’s, the load on the wind tunnel model must 
be scaled to be the value W„ using the relationship: 


Wyy 

^ref,w 



CLqooj 


(125) 


However, the wing loading on the wind tunnel model must be adjusted to take into account the change in dynamic 
pressure to preserve the lift coefficient Cl- Multiplying Lq. 125 by the dynamic pressure ratio j . a new value for 
the load for a design Cl can be determined. 


^ ^ ^ = CLq^J ( ^ ) 

^ref .w ^ref \ / \ / 


Cl^oo^w 


(126) 


With respect to the torsional stiffness, the analysis beginning with Lq. 100 is still valid. However, the assumption 
that the dynamic pressure value qoo remains constant in Lq. 106 is no longer valid. Instead, the analysis needs to be 
adjusted as follows: 


“ {cm^Cf + CLefCOsA)qoojCOS^ASeJ _ - {c„,^^c„ + ClCwCOsA) q^^„cos^ ASe^ 


GJf 


GJ„ 


(127) 


f clc f cos A'j q^o^f cos AS^ j (c^ac 3 ~cl cos ^oo.wcos AS^jjn^ 


2 

scale 


GJf GJ„ 

- {cmacC f + CLef cos A) qoojcos^ASej _ - {cm^c f + CLCf cos a) qoo, fCOS^ASeJ f q. 


GJ 


f 


^■^'vn^scale 


oo,W 


GJ„ = — GJf 

'‘'scale 

The bending analysis is also altered from Lq. 122. 


Jooj) 


(128) 

(129) 

(130) 


ClJooJ COS^ ASejH fY CLqoo,„ COS^ ASe,w{lwY 

'■“» —f 

CLqoojCOS^ ASej{lfY _ CLjoc, wCOS^ A{SeJ (If/nscalef 
bending rFr~\ ~ TFj~\ 

\^^yy)w 


^bending 


{EIyy)f 

ClCOS^ ASej{lfY ClCOS^ ASgj{lfY / qoo. 


{EIyy)f 


{EIyy)wtl scale 


qooj 


(131) 

(132) 

(133) 
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1 


(134) 


{EIyy)w — 


nbendingn^^alg 


-{EIyy)f 




It can be seen that in test situations where the wind tunnel dynamic pressure cannot be the same 
flight condition, an additional factor equal to the dynamic pressure ratio j must be added to the 
factors r,orsion and rbending such that 


1 

^ torsion — 

11'^ 

'‘'scale 




as the full-scale 
stiffness scaling 

(135) 


1 

^bending — 4 

nbendingn scale 


/ qoo,w 


(136) 


4. Mach Number Ejfects 


The above analysis considers the case where the change in the dynamic pressure from qooj to qoo,w does not affect the 
aerodynamic coefficients. That is, = ^mac and ctj = c/,,w = cl- This assumption allows for the clean 

derivation of the analyses above. In reality, the change in the dynamic pressure affects the Mach number of the flight 
condition, which affects the aerodynamic coefficients. 

If the wind tunnel is operating at a different dynamic pressure q^o^w and at ambient sea level altitude, the Mach 
number of the wind tunnel model can be calculated. 


My^ 


/ 2^00, W \ 2 

V ypsL ) 


(137) 


This Mach number is different than the Mach number Mf of the full scale wing alone model and is expected 
to be lower. This affects the aerodynamic coefficients. 

The aerodynamic data for the wind tunnel model is obtained for two different flight conditions, assuming a rigid 
undeformed wing planform. 

• The ESAC’s cruise condition at Mf = 0.797, h = 36,000 ft. For Cl = 0.510 the angle of attack is taken to be 
a = 2.249°. 

• An estimated wind tunnel flight condition of qooj = 20^, sea-level flight, corresponding to a Mach number 
M„ = 0.1 16. For Cl = 0.510 the angle of attack is taken to be a = 3.993°. 


The lift curve of the wind tunnel model at both flight conditions are compared. The results show that the total aircraft 
Cl„ is different as a result of the Mach number effect, which is to be expected due to compressibility. A higher angle 
of attack for level flight is required to achieve the same Cl = 0.510 at the lower Mach number as well. 



Figure 15. Mach Number Effect on Lift Curve 
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The spanwise lift c/ and moment c„i^^ distributions are shown in Fig. 16. 




Figure 16. Mach Number Effect on Spanwise Aerodynamic Coefficients 


The spanwise lift distributions are similar in magnitude due to the fact that the overall load is maintained at Cl = 
0.510, but it is observed that the load shifts slightly outboard at the higher Mach number. The spanwise distribution 
of c,n^^ is drastically different between the flight conditions. At the lower Mach number, the values of are much 
lower. Thus, the aeroelastic scaling for the wind tunnel model will require consideration that Mach number effects can 
alter the spanwise Cmac ■ 

The scaling analysis conducted for bending still holds because does not affect the derivation of Eq. 136. 
However, additional analysis is required for aeroelastic scaling of torsion. The analysis is modified from Eq. 127. 
While the simplification CL.f = cl^w = cl is still made, Cmacj ^ ^ntacw- 


- {cmacjCf + Cie f cos A) cos^ AS^j 

GJf 


“fCr. e? videos A) AS^^ve 


(138) 


“ + CLBf cos A) qoojcos'^ ASeJ _ - {cmacW [cf/n^cale] +CL [ef/riscale] COS A) q^o^wCOS^ A {Se,w/n^cale) 


GJf 


GJ„ 


- + Ctef cos A) qoojCOS^ ASej _ - {cmacwCf + CLef cos a) qoojCOS^ ASeJ ( qoo^w 


GJf 




GJ„ = GJ 


„3 


1 


scale/ VCmr./ 


Zac I 

qoefcos ASefCfl 

q«,j 


where 


and 


Cmr.f = Cm^cJ + COS A = ( Cm^c + COS A 
C/ V Cf 


Cmfw — w ' 


Cf 


cos A = c. 


cief 

Cf 


cos A 


(139) 

(140) 

(141) 

(142) 

(143) 


Eor the flight condition, c„, is generally a negative number and is generally a positive number. However, the 
Mach number effect introduces the possibility that c„r,f and c,„rM are opposing in sign due to the relative magnitudes 
of the terms. This causes the ratio < 0- In effect, reducing the Mach number from Mf to M„ can cause the 

wing aerodynamics to attempt to twist in an opposite direction. If this results, then GJ„ becomes a negative number, 
which makes no physical sense. The interpretation is that there exists value of M„ or wind tunnel dynamic pressure 
qoo^w at which the ratio j is negative and scaling the torsional stiffness of the wind tunnel model cannot achieve 
the same amount of twist as the full scale flight condition. 
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The values c,„r,f and c„,r.w are related to the torsional force experienced by a wing element about the elastic axis. 
For the ESAC’s cruise condition at Mf = 0.797 and the wind tunnel model’s Mach number for qocj = 20^, sea-level, 
the values of c„,r are compared. 



Figure 17. Mach Number Effect on Torsional Force c„r At Q, = 0.510 

It can be seen that Cmr.f and c,nr,w are the same in sign for most locations along the span. For the purpose of scaling 
torsional stiffness, this allows it to be possible for the wind tunnel model’s tip twist to be scaled to reflect the full-scale 
wing alone model’s at cruise. 

For the purpose of illustration, the values of c,nr for two other flight conditions are plotted in Fig. 18. The two 
flight conditions correspond to a C/, = 0.346. The full scale wing alone corresponds to Mf = 0.8 at a h = 30,000 ft 
altitude. The wind tunnel flight condition corresponds to qocj = 20^ at sea-level flight corresponding to M„ = 0.1 16. 



Figure 18. Mach Number Effect on Torsional Force c^r At Q, = 0.346 

It can be seen that Cmrj and Cmr,w are opposite in sign for all locations along the span when scaling from Mf = 0.8. 
This indicates that the ratio is generally negative along the span. The Mach number effects prevent the wind 
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tunnel model’s tip twist at the wind tunnel condition to be scaled to the cruise M = 0.8 condition. 

5. Coupled Model Considerations 

The above analysis is intended for models that do not possess a coupled aerodynamic-structural nature and can be 
readily applied to problems where bending and torsion do not exhibit strong coupling. The actual static aeroelastic 
problem for the wind tunnel model possesses both coupling with an aerodynamic model, and the structural bending 
and torsion modes are also coupled together as a result of aerodynamic stiffening and softening. 

Development of an analytical estimate for a scaling factor for the coupled model can be extremely intensive. 
Instead, additional factors ri and r 2 representative of scaling factors are added into the relationships that should be 
tuned heuristically using the static aeroelastic framework in Fig. 13 for the design Cl = 0.510. For stiff wing problems, 
the coupled results and the uncoupled results are generally close and the scaling factors determined for the uncoupled 
model can be used. As even softer wing models are developed with higher vertical tip deflections, the usage of ri and 
T 2 factors become more necessary. 


GJ„ = 


^torsionGJ f — r\ 


(^) GJf 

Kale 

- 

f torsion 


(144) 


{EIyy)„ = rbending{EIyy) f = T2 ^ ) (Elyy) f (145) 

n bend itign scale V / 

' V ^ 

^bending 

The factor r\ is also used to represent a constant scaling factor to correct for mach number effects in Eq. 141. 

The relative percentage of aeroelastic wing tip deflection and the relative aeroelastic twist about the elastic axis 

20 

are determined using the static aeroelastic model of scaled wind tunnel model using the scaling factors rtorsion 
and rbending in Eqs. 144 and 145. The results are summarized in Table 5 and 6. 


100(^),% 

T2 = 0.50 

t2 = 0.80 

T2= 1.00 

r 2 = 1.50 

ri = 1.00 

14.474 

10.032 

8.343 

5.871 

n = 1.50 

14.465 

10.025 

8.335 

5.868 

n = 1.75 

14.466 

10.020 

8.336 

5.868 

ri = 2.00 

14.468 

10.020 

8.333 

5.867 


Table 5. Aeroelastic Relative Wing Tip Deflection for Wind Ttinnei Model with ri and V 2 Scaling Factors 


/ V deg 

^ IT 

t2 = 0.50 

r2 = 0.80 

O 

p 

II 

t2= 1.50 

ri = 1.00 

-5.577 X 10-3 

-7.251 X 10-3 

-7.923 X 10-3 

-8.958 X 10-3 

O 

II 

C 

-3.717 X 10-3 

-4.828 X 10-3 

-5.264 X 10-3 

-5.965 X 10-3 

II 

-3.195 X 10-3 

-4.129 X 10-3 

-4.515 X 10-3 

-5.113 X 10-3 

n = 2.00 

-2.789 X 10-3 

-3.613 X 10-3 

-3.946 X 10-3 

-4.472 X 10-3 


Table 6. Aeroelastic Relative Wing Tip Elastic Axis Twist for Wind Ttinnei Model with r\ and T 2 Scaling Factors 


6. Summary of Aeroelastic Scaling 

The scaling equations represented by Eqs. 144 and 145 are applied to the stiffness distributions of GJ and Elyy. While 
rtorsion and rbending Can be defined as functions of span, this would result in a more complicated scaling procedure. 
Instead, the values of rtorsion and rbending are selected as single values applied to the entire baseline GJ and Elyy 
distributions. This simplification can be done because the wind tunnel model’s aeroelastic behavior is scaled based on 
its tip deflection and tip twist. 
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Given a value of riscaie = 10.3585, aeroelastic scaling of the wind tunnel model from the cruise condition of Mach 
number Mf — 0.797, h = 36,000 ft down to qoo,w = 20^ (M„ = 0.1 16), sea-level is conducted, where the design lift 
coefficient is = 0.510. The scaling factors that were determined are summarized in Table 7. 


Cl 

0.510 

^bending 

0.537 

n 

1.75 

ri 

0.80 


9.47 X 10-2 

^bending 

1.28 X 10-"^ 

^ torsion 

4.42 X 10-® 


Table 7. Aeroelastic Scaling Factors For Wind Tnnnel Model 


The application of the determined scaling factors results in aeroelastic deflections for the flexible wind tunnel 
model summarized in Table 8. The scaling factors are able to scale the wind tunnel model such that the wing tip 

deflection is 10.02% relative to the wind tunnel model’s semi-span, and the relative elastic axis twist is ) = 


-4.129 X 10-3^, 


which compares to the full-scale model’s elastic axis twist of 


(^)^ = -4-135 X 10-3^51. 


Cl = 0.510 

Lull Scale Wing Alone 

Wind Tunnel Model 

a, deg 

3.064 

5.773 

W,ip, ft 

2.740 

0.540 

lOO(^), % 

4.907 

10.020 

®tip, deg 

-0.231 

-0.223 X 10-1 

/ ^®tip \ deg 

-4.135 X 10-2 

-4.129 X 10-2 


Table 8. Aeroelastic Deflection Results for Wind Tunnel Model with Aeroelastic Stiffness Scaling 


C. Static Divergence 

An analysis of the scaled down wind tunnel model is conducted to determine the divergence dynamic pressure qd, or 
the dynamic pressure in which the wind tunnel model will experience static divergence. Determining the divergence 
dynamic pressure places a restriction on the wind tunnel test condition and is important to analyze to ensure that the 
model will be able to be properly utilized at the wind tunnel test conditions. If the divergence dynamic pressure is 
significantly larger than the test condition of the tunnel, static instability does not pose a problem. This, however, does 
not preclude the possibility of the wind tunnel model experiencing dynamic instability due to aeroelasticity, or flutter. 
Flutter is not investigated within the scope of this study. 

1. Torsional Divergence 

Initially, a preliminary analysis of static divergence can be performed on the wind tunnel model focusing only on 
torsional divergence. Torsional divergence is a classically examined phenomenon due to the basic aeroelastic coupling 
in twist.^’ Aerodynamic forces can cause a wing to twist nose-up (negative 0). This nose-up twist increases the 
aeroelastic angle of attack on the wing sections. Since lift force is proportional to aeroelastic angle of attack, this 
can cause the wing to twist even more nose-up. Thus, a positive feedback loop exists between twist and angle of attack 
that can cause the wing to exceed its structural limitations and experience torsional divergence. 

The structural stiffness and aerodynamic stiffness matrices for the globally assembled finite-element system were 
previous represented as Ka and in Eq. 85 and assembled by the element matrices in Eqs. 80 and 81. Let the global 
stiffness matrices be partitioned as follows: 



^s,tb 

^sM 

Ks^b 
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(147) 


K„ = 


Ka,t 

^a,tb 

Ka.bt 



where Ks^ and A), ; are sub-matrices of the stiffness matrix elements corresponding to the torsional degrees-of-freedom 


01 02 


0fi+i 


, and A), j, are the sub-matrices of the stiffness matrix elements corresponding to the bend- 


ing degrees-of-freedom 


Wl Wj W2 


Wn+\ 


''n+\ 


. and Ksjb, K,h„ Katb, and Kaht are 


the coupling matrices. 

A torsional divergence analysis involves examining the matrices Aj; and A),;. Four cases are examined using 
different scaling factors rtorsion- 


• A first scaling case where rtonion = ri - 4 ^ ^ = 1 .00, n^caie = 10.3585, and = 9.47 x 10 

"scale V?”.// ■ 

• A second scaling case where rtorsion = r\ ,-^ = \ qo, tiscaie = 10.3585, and = 9.47 x 10^^. 

"scale ^ 

• A third scaling case where rtorsion = ri , ri = 1 .75, tiscaie = 10.3585, and =9.47 x 10^^. 

"scale ^ 

• A fourth scaling case where rtorsion = ri f ^ j riscaie = 10.3585, and =9.47 x 10^^. 

"scale ^ 


The aerodynamic stiffness matrix Ka^t is reliant on dynamic pressure as seen in Eq. 81, and static divergence occurs 
when the term Ai, , -f A), ^ is non-invertible or singular. In order to evaluate this, let the determinant of the total stiffness 
matrix normalized to determinant of the zero-speed structural stiffness matrix be defined as 


det(/f5 , -f Kot ) 
det(/:i,,) 


(148) 


The wind tunnel test facility is limited at operating dynamic pressure q^o < 60^. For this analysis, however, 

the value of A, is plotted versus ranging up to qoe = 400 for illustration purposes. The value in which A, = 0 
represents the divergence speed qi}. 



Figure 19. A, versus Dynamic Pressure q^. Torsional Divergence Analysis 


A torsional divergence dynamic pressure of ^^ = 162^ is observed for the first scaling case where rtorsion is the 
smallest of the four scaling cases. The third scaling case has a torsional divergence dynamic pressure of q^ = 274|^, 
while the second and fourth scaling cases where rtorsion have divergence dynamic pressures beyond q„ = 

^^xcale 
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400^. It is clear that increasing the value of Vtorsion, and thus the torsional stiffness of the model, causes the torsional 
divergence dynamic pressure to increase due to the fact that the additional structural stiffness helps to prevent onset 
of structural instability. For all four scaling cases, the torsional divergence dynamic pressure is far beyond the desired 
wind tunnel test speed of = 20|^, and thus, torsional divergence is not expected to be a problem for the test 
condition of the sub-scale wind tunnel model. 


2. Coupled Static Divergence 


An analysis of torsional divergence represents an uncoupled analysis of the static instability problem for aeroelasticity. 
The full FEM model developed in this study can actually be used to conduct a more refined analysis by examining 
the full Ks and matrices in Eq. 85. Based on the sign of the terms in Eq. 81, aeroelasticity is expected to result in 
softening in torsion and stiffening in bending. The bending slope contributes to the aeroelastic angle attack through 
Eq. 46 and relieves the angle of attack, actually improving the divergence properties of the aeroelastic model and 
indicating that sole analysis of torsional divergence can actually be more conservative than the real system. Static 
divergence analyses for four scaling cases are examined: 

• A first scaling case where r,orsion = n T— ( ) > fbending = 4? — ( ) ’ '‘i = ^ '‘2 = 1 -00^ 'Aea/e = 

"scale J * "bending" scale J 

10.3585, and ^ = 9.47 x 10^1 

• A second scaling case where r,orsion = n f ) . fbending = 4? — f ) > ''i = 1 -00, r 2 = 1 .00, ngcaU = 

^ scale s / "bending^scale \ / 

10.3585, and ^ = 9.47 x lO^^. 

• A third scaling case where rtorsion = n ) > lending = ^— 4 — ( ) > ''1 = 1-75, ^2 = 0.80, = 

10.3585, and ^ = 9.47 x lO^^. 

• A fourth scaling case where rtorsion = n f ) . ''feenrfmg = '— 4 — ( ) , ri = 1 .75, r 2 = 0.80, ngade = 

'^scale s / "bending^scale s / 

10.3585, and ^ = 9.47 x lO^^. 

9 -,/ 


Let the determinant of the total stiffness matrix normalized to the determinant of the zero-speed structural stiffness 
matrix be defined as 


det(jg;,-fjr„) 

det(fi'^) 


(149) 


The value of A is plotted versus qoo to determine the value qd when Kg + Ka becomes singular. 



Figure 20. A versus Dynamic Pressure q^. Coupled Static Divergence Analysis 
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The results in Fig. 20 demonstrate the effect of modeling coupled bending-torsion on static divergence. The term 
A is positive for all the scaling cases within the dynamic pressure range examined. In fact, addition of the coupled 
bending-torsion consideration actually causes the term A to increase as dynamic pressure increases, indicating that 
static instability is not an issue for the model. 

Let the total stiffness matrix be represented as 


K = K, + Ka 


(150) 


This can be expanded using Eqs. 146 and 147 as 


K = 

Kt 

Ktb 


Kg,, 

^s,tb 

+ 

Ka,t 

^a,tb 


_ Kbt 

Kb 


_ Kg,b, 

Ks,b _ 


Ka,bt 

^a,b 


(151) 


The terms K^^tb and K^^bt are negligible because B\ and B 2 in Eq. 18 are considered to be zero. Thus Eq. 151 
becomes 


K, 

K,b 


Kg,, + Ka,, 

1 

_ Kb, 

Kb 


1 

Sr 

^s,b “1" ^a,b 


(152) 


The matrix in Eq. 152 can be expressed as 


Ks,t + Kaj Ka^tb 

Ka^bt Kg b + K(i b 


I Ka^tb 
0 Kgb + Kafi 


Kg,t + Ka, - Ka^tb {Kg,b + Ka,b) ^ K, 


[Kg,b + K,^b) ^Ka^bt 


a.bt 


0 

I 


(153) 


which allows the determinant to be calculated as 


det(/r) 


det(A"^^i -f A"aji)det(/T^ ( -f A), f A), -f A), K^ bt) 

det( A',,^ + Ka^b)A&t{Kg^t + Ka,t - Ka^tbR-lK-^bt) 


(154) 


Because bending stiffness does not experience static instability, the determinant det(A^, -f A^, £,) is positive. While 
it is difficult to make any generalizations about the sign of the second determinant term in Eq. 154, it can be concluded, 
however, that static divergence occurs only if 

det(A',., + Ka.t - Ka,tbK^},K-lKa,bt) <0 (155) 


It is possible that the hnal term in Eq. 155 can be always positive, det(A^ , -f Ka,t — Ka^tbK^l^^bKafit) > 0, based 
on the values of the elements in the stiffness matrices. This is the case for the aeroelastically scaled wind tunnel model 
whose results are in Eig. 20 where static divergence does not occur for the model. It can be seen that the torsional 
divergence problem is alleviated based on the —Ka^tbK^l^^bKa,bt term in Eq. 155. 


IX. Wing Twist Optimization 

To complete the development of the wind tunnel model, the geometrically and aeroelastically scaled model needs 
to be re-twisted for the wind tunnel test condition. A new unloaded shape is developed such that when the flexible 
wing model is operating at wind tunnel test condition, it aeroelastically deforms to a deflected shape that has minimum 
induced drag or maximum L/D ratio. This tailors the model such that it becomes ideal for conducting trade studies 
and drag analysis. An optimization procedure is developed and applied to the sub-scale wind tunnel model. 


A. Optimization Method 

Optimization is achieved using an unconstrained gradient-based optimization algorithm. The foundation of a gradient- 
based optimization method is the determination of an optimal search direction which sufficiently minimizes the objec- 
tive function, calculated using the function gradient information. Eor this particular problem, the objective function 
is not an explicit analytical function, but instead the static aeroelastic mapping developed in Eig. 13. The input flight 
condition for the wind tunnel model is hxed, where Cl = 0.5 10, Mach number M = M„ = 0.116, and altitude h = 0 
corresponds to sea-level testing conditions. However, a new design input is added which allows a user to add addi- 
tional twist onto the sub-scale model’s existing pre-twist distribution. Let the design input be expressed as a x, which 
contains t] individual variables that specify an additional twist distribution on the sub-scale model. 
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The models utilizes Xj to control the wing twist shape, and the static aeroelastic model determines the aeroelastic 
shape of the (©, W, Wx) and the angle of attack a required for the static aeroelastic model to have a lift coefficient 
equivalent to Cl- Aerodynamic modeling of the deformed shape also allows for the drag coefficient of the model Co, 
to be determined for the static aeroelastic model. For all intents and purposes, then, the static aeroelastic model can be 
seen as an equivalent functional mapping such that 




(156) 


where the objective function Jc is accomplished by utilizing the static aeroelastic mapping. 

Because Jc is not expressed analytically, the gradient of the objective function cannot be explicitly calculated and 


must be approximated. Let Xi = 


. A forward hnite-difference method is used to approximate the 


gradient about a known design point, x, p using the following: 


VJc{Xi) 

dJcjxifl) 

dxifl 


dJc dJc dJc 

dxj dxf • ■ ■ dx^ 

Jc(Xifi + AXi)-Jc{Xifl) 

Axi 
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The gradient of a function points in the direction of greatest increase, so a possible search direction is one in the 
exact opposite direction of the gradient itself. This is called the direction of steepest descent. However, while choosing 
steepest descent will result in convergence on a minimum, it is known to be slow and inefficient.*® Therefore, several 
other methods have been developed to determine a more efficient search direction, such as the method of conjugate 
directions.*^ The conjugate direction method uses the following property of conjugate vectors 


Si[H]Sj = 0,i^j (159) 

where [H] is a symmetric and positive-dehnite matrix. For a quadratic function, where [//] is the Hessian of the 
function, if conjugate search directions S are used, then the conjugate property produces a complete decoupling which 
results in the ability to optimize the function in exactly rj line searches, where Tj is the number of problem design 
variables. While this method is the most efficient for exactly quadratic problems, it is also effective for non-quadratic 
functions. 

To start the conjugate direction method, the first search direction, 5**, is calculated using steepest descent since there 
is no prior gradient information. For following iterations, each updated search direction is found using information 
about the gradient at the current design point and the search direction from the previous iteration. The formula for 
updating the search direction at each iteration is 


= -VJ^ + /5S^-^ 

where 


(160) 


(V7,^)^(V7,^-V7^*) 

(V7^*)^(V7^*) 


(161) 


using the Polak-Ribiere method,*^ one of several different possible j3 formulations. 

It is possible for the search direction to become ill-conditioned due to numerical imprecision or if the objective 
function is particularly non-quadratic. Therefore, it is necessary to monitor whether or not the conjugate search 
direction does not result in a sufficient decrease in the objective function and reset to steepest descent if necessary. 
There are two particular scenarios in which the search direction should be reset. The hrst situation occurs when the line 
search does not produce an improvement in the objective function. Second, each time a conjugate search direction is 
calculated, it should be compared with the gradient of the function at that design point. The closer the search direction 
is to the direction of the gradient, the less likely it is to result in sufficient decrease. In particular, the angle between 
the conjugate search direction vector and the negative of the gradient should be less than 90°. 

With the search direction determined, a one-dimensional optimization in the search direction, or line search, is 
conducted. For the line search, the design variable is the search step in the particular search direction. The formal 
statement of the problem is as follows: 


min7c(jcf + t > 0 

t 
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where jcf are the actual problem design variables at iteration step k, is the search direction at iteration step k, and <7 
is the search step variable for the line search. 

The line search consists of several steps in order to determine a minimum of the objective function in the search 
direction. First, basic bracketing*® is used to determine a lower and upper bound between which a minimum in the 
objective function exists. With upper and lower bounds defined, the bracketed interval is further refined using the 
Golden-Ratio Search method.*® This method uses the Golden-Section ratio value to reduce the interval around the 
minimum to a desired tolerance, in the fewest number of function evaluations. 

For this method, two interior points in the interval are calculated using the following formulas: 

(Ja — (yiwr~^ '^{^upr ^Iwr) ( 163 ) 

Ob = (J/w + ( 1 - T:){Oupr - Oi„r) (164) 

where T is derived from the Golden-Section ratio, and is given by 


t = 3 ^ = 0.3819 
2 

The objective function is then calculated at the two additional interior points, so that all four points. 


(165) 


C/m 


Oa Ob Oi 


upr 


, and their corresponding function values. 


J, 


cjwr 


Jc.a 


c,b 


J, 


c,upr 


, are known. The 


interval refinement algorithm then proceeds as follows: 


• If Jc,b > Jc.a, then Ob becomes the new upper bound, Oupr, and Oa becomes the new o* interior point. A new 
interior point, Oa, is calculated using Eq. 163, along with it’s corresponding function evaluation, Jc.a- 

• If Jc.b < Jc.a, then (Jfl becomes the new lower bound, Oi^r, and o* becomes the new Oa interior point. A new 
interior point, Ob, is calculated using Eq. 164, along with it’s corresponding function evaluation, Jc.b- 

This process is continued until the interval has been reduced to a desired level of accuracy relative to the original 
interval. 

Once the final interval bracketing the minimum has been calculated, the minimum of the objective function in 
the particular search direction can be approximated using a polynomial fit. In this case, since there are four points 
available from the Golden-Search method with four known function values, a cubic polynomial approximation of the 
objective function can be determined from which the approximate minimum can be calculated by finding the roots of 
the derivative of the polynomial approximation. 

The result of the line search is the minimum of the objective function in the particular search direction and the 
corresponding minimum search point, <7*, which is then used to find the next design point as follows: 


xf+* =xf-b (7*5^ (166) 

This design point is then used as the initial design point for the next iteration of the optimization. 

The method of calculating the search direction, performing a line search, and updating the search direction at the 
new design point is continued until a convergence criteria is met. Eor this problem, convergence is assumed when the 
absolute value of the objective function has not changed over several iterations. 


B. Design Variable Distributions 

The design variables for this particular wing twist optimization problem control the values of the additional twist to be 
added to the already existing jig-shape twist of the flexible wing. Let the wing pre-twist y{x) be the existing pre-twist 
on the jig-shape of sub-scale wind tunnel model, positive nose-down, about the pitch axis of the wind tunnel model. 
Let the total wing pre-twist y be represented as 

y{x) = 7 (x) -f A 7 (x) (167) 

where A 7 (x) represents an additional pre-twist, positive nose-down, applied to the jig-shape on the wing about the 
pitch axis. The design variable x, controls the distribution of A 7 such that 

Ay = f(xi) (168) 

Two different design variable cases are considered for the wing twist optimization: discrete point, and polynomial 
shape function coefficients. 
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1. Discrete-Point Design Variables 


For the first case, the design variables to be input into the optimization method are the actual additional twist values 
at two discrete points on the wing. In particular, the points are at the wing break point due to the wing trailing edge 
extension x^reak, and the wing tip xnp. That is, for a discrete-point design variable optimization 


Xi = 


^Yi^break) ^Yi^fip) 


(169) 


The additional twist is also constrained to be zero at the wing root A7(0) = 0. The total twist distribution y{x) is 
determined by linearly interpolating A 7 (x) across the wing span and adding it to y{x). 


2. Shape Function Design Variables 

For the second case, the additional twist along the wing span is represented by Chebyshev polynomial functions. In 
particular, the following shape function is initially considered; 


A 7 (x) = ao7b(x) -|-ai7'i(x) -|-a 2 T 2 (x) -|-a 373 (x) -\-a/\T 4 {x) (170) 

where 

7b = 1 (171) 

Ty=X (172) 

72 = 2x^-1 (173) 

73=4x^-3x (174) 

74 = Sx^* - 8x2 -7 1 


However, in order to compare directly with the discrete optimization and prevent a under-constrained optimization 
problem, the functions needs to be modified such that the additional twist is always fixed to be zero at the wing root, 
A7(0) = 0. This is done by subtracting the root value from Equation (170). For this particular model, the constant 
terms in the equation are eliminated. 

Additional scaling of the shape function is done so that the polynomial coefficients stay nearby in order of magni- 
tude. The location along the wing x is scaled by the length of the wing L. Therefore, the final shape function that is 
used to describe the wing twist distribution for the model is given as 


The design variables for the shape function optimization are the four polynomial coefficients of the above shape 
function 

a\ aj a^ a^ (177) 


C. Optimization Results 

The optimization framework is utilized to determine the new wind tunnel undeformed shape, such that when aeroelas- 
tically deformed at Cl = 0.510, Cd is minimized. This corresponds also to an optimized L/D for the wind tunnel test 
condition. 

7. Discrete-Point Optimization 

A series of optimization runs are conducted to minimize Cd using wing pre-twist specified at discrete points along the 
wing. The value of the additional pre-twist is fixed such that no additional wash-out is added to the root A7(0) = 0, 
but the pre-twist at two locations on the wing are prescribed as design variables and linearly interpolated for stations 
in between. The two locations selected are the wing tip located at = 5 .388 ft, and y break = 1-793 ft along pitch axis 
or the 7 b— axis in Fig. 12 or the l? 2 ~direction in Fig. 8 . 

Without adding any additional pre-twist, the wind tunnel model has aCo = 0.03228. Four optimization runs are 
conducted using the discrete-point design input, each initialized at different starting values. The total amount of design 
variables for the discrete-point optimization is Tj = 2 corresponding to A 7 at the wing extension break and the wing 
tip. 
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The results of the value of the cost function Jc or Co of the four optimization runs which were initialized at different 
starting values, are plotted in Fig. 21. For this case with only rj = 2, the optimization requires only a minimal amount 
of iterations before converging to the minimum Cjj value. 



Figure 21. Cost Function (Cn) per Iteration of Discrete-Point Optimization 

The optimization result corresponded to Cd = 0.03020, representing a 6.444% decrease in Co relative to the un- 
optimized wind tunnel model. The additional pre-twist distribution that is added to the wind tunnel model planform 
results are summarized in Table 9, and a plot of the pre-twist distribution along the wing shown in Fig. 22. It can be 
seen that a nose-up additional pre-twist Ay is imposed on the wind tunnel model in order to reduce the Cd value at the 
wind tunnel test condition. 


y, BBL, ft 

Ay, deg, positive nose-down 

0 

0 

1.793 

-4.220 

5.388 

-6.203 


Table 9. Optimization Result for Discrete-Point Optimization 



3 

y. BBL, ft 


Figure 22. Optimized A/ Result for Discrete-Point Optimization 
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Deflection information for the static aeroelastic model for the un-optimized and the discrete-point optimized results 
are summarized in Table 10. 


Cl = 0.510 

Un-optimized Wind Tunnel Model 

Optimized Wind Tunnel Model 

a, deg 

5.773 

2.411 

Cd, counts 

322.8 

302.0 

%p, ft 

0.540 

0.649 

lOO(^), % 

10.020 

11.872 

@tip, deg 

-0.223 X 10-1 

-0.343 X 10-1 

/ ^®tip \ deg 

^~b~ ^ IT 

-4.129 X 10-3 

-6.371 X 10-3 


Table 10. Aeroelastic Deflection Results for Discrete-Point Optimized Wind Tlinnel Model 

The lift distribution of the un-optimized and discrete-point optimized wind tunnel models undergoing aeroelastic 
deformation are also shown in Fig. 23. 



Figure 23. Lift Distribution of Wind Tunnel Model Using Discrete-Point Optimization Results 

The un-optimized model’s lift distribution is generally triangular in shape, due in part to the aeroelastic deformation 
which causes the wing tip to twist nose-down. This nose-down aeroelastic deformation causes the lift distribution to 
shift towards the wing root. The addition of the optimized Ay corrects this, however, and a new pre-twist is prescribed 
that twists the wing tip more nose-up as shown in Fig. 22. The resulting lift distribution becomes more elliptical in 
shape. 

2. Shape Function Optimization 

A second series of optimization runs are conducted using the shape function in Eq. 176 to prescribe additional pre- 
twist on the wing Ay. In this case, the design input variable x, has increased to rj = 4 degrees -of-freedom, where the 
input design variable represented in Eq. 177, corresponds to the shape function coefficients in Eq. 176. The shape 
function imposes no additional wash-out added to the root, or A7(0) = 0. 

A total of seven optimization runs are conducted from different initial values, and the evolution and decrease in 
the value of Cd as iterations are conducted is plotted in Eig. 24. While it requires more iterations before the cost 
function decreases to a minimum value, Eig. 24 shows that the optimization is able to drive Cd of the wind tunnel 
model at the test conditions down. The independent optimization runs also generally converge to similar minimum 
Cd values. The lowest minimum Cd value obtained in the optimization study was Cd = 0.03018 corresponding to a 
6.506% improvement from the un-optimized wind tunnel model shape. The shape function optimization Cd result is 
lower than that of the discrete-point optimization, but is to be expected due to the fact that Ay is parametrized by more 
degrees-of-freedom in the shape function optimization. 
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Figure 24. Cost Function (Co) per Iteration of Shape Fnnction Optimization 

The optimization result is summarized in Table 1 1 for the parameters of the shape function, where it is seen that 
normalization of the independent variable in Eq. 176 is able to return coefficient results of similar order. 


Parameter 

Value 

ai 

-0.3379 

fl2 

-0.9700 


0.5701 

G4 

1.3101 


Table 11. Optimization Result for Shape Function Optimization 

The optimization results translate into an additional pre-twist distribution Ay{x) added to the wind tunnel model 
planform, and the distribution is shown in Fig. 25. By using polynomial basis functions, the pre-twist distribution is a 
smooth continuous curve. 



3 

y. BBL, ft 


Figure 25. Optimized Ay Result for Shape Function Optimization 

The static aeroelastic model is used to obtain the deflection results for the un-optimized and the shape function 
optimized results, and the values are summarized in Table 12. 
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Cl = 0.510 

Un-optimized Wind Tunnel Model 

Optimized Wind Tunnel Model 

a, deg 

5.773 

2.345 

Cd, counts 

322.8 

301.8 

Wnp, ft 

0.540 

0.644 

100(2^),% 

10.020 

11.949 

@tip, deg 

-0.223 X 10-1 

-0.355 X 10-1 

/ ^®tip \ deg 

\~IT~ J’ “FT 

-4.129 X 10-2 

-6.594 X 10-2 


Table 12. Aeroelastic Deflection Results for Shape Function Optimized Wind Ttinnel Mode 

The lift distribution of the un-optimized and shape function optimized wind tunnel models are also shown in Fig. 
26. The lift distribution of the wind tunnel model with wing pre-twist optimized using shape functions is very similar 
to the lift distribution of the model with wing pre-twist optimized using the discrete-points. Both optimization results 
impose a nose-up twist onto the wing to counteract the nose-down twist due to aeroelastic deformation. 



Figure 26. Lift Distribution of Wind Tlmnel Model Using Shape Function Optimization Results 


3. Summary 

The results of the discrete-point optimization and the shape function optimization are used to analyze the static aeroe- 
lastic model of the wind tunnel model. Lift curves and drag polars are plotted in Fig. 27 representing the flexible wind 
tunnel model, and the curves provide insight into the optimization results of the model. 

Both the discrete-point optimization result and the shape function optimization result produce lift curves which 
are very similar to each other, and the lift curves for the optimized models are shifted upwards of the un-optimized 
lift curve. This means that the additional optimized pre-twist Ay helps to recover the loss of lift due to the nose-down 
aeroelastic deformation of the flexible wind tunnel model. 

The drag polar of the flexible sub-scale wind tunnel model is also affected when the optimized pre-twist Ay is 
applied to the static aeroelastic model. Though slight, the drag polars of the optimized models are shifted, and the drag 
polar at the design Cl value is lower than that of the unoptimized model. To recall, the discrete-point optimized model 
observed a 6.444% decrease in Cd at the test condition relative to the un-optimized model and the shape function 
optimized model observed a 6.506% decrease. 
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Figure 27. Lift Curves and Drag Polars for Aeroelastic Wind Ttinnel Models 
The aeroelastic deflections for the optimized and un-optimized aeroelastic models are plotted in Fig. 28. 




Figure 28. Aeroelastic Deflections for Wind Ttinnel Models 

Because of the nose-up pre-twist applied that resulted in the shift of the lift curves of the optimized models, the 
vertical bending deflection Wtip for the optimized models is higher. The higher lift also contributes to the moment 
about the elastic axis of the wing, thus driving ©„p more negative, or nose-up. 

The aircraft rigid body angle of attack is a, the aeroelastic deformation effect on the angle of attack is Ue, and Ay 
represents the additional prescribed wash-out determined through optimization-all about the aircraft pitch axis where 
a and Ue are positive nose-up, and Ay is positive nose-down. Let a new quantity Up be defined such that 

ap = a + ae-Ay (178) 

where Up represents the angle of attack of a local section of the wing perpendicular to the pitch axis, positive nose-up, 
relative to the wind tunnel un-optimized jig-shape existing pre-twist y. Thus, Up represents the effective angle of attack 
of a wing section relative to the un-optimized existing pre-twist y. The value of (Xp^tip is plotted versus the aircraft 
angle of attack a for the aeroelastic wind tunnel models in Fig. 29. 
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Figure 29. Physical Angle of Attack at Wing Tip for Aeroelastic Wind Ttinnel Models 

The plot of Ctp shows the physical angle of attack of the wing tip relative to the jig-shape for the un-optimized 
and optimized models, where it is the largest for the model optimized using the shape function. For an un-optimized 
model with a rigid planform, «p = a, and thus the plot of dp, tip for the un-optimized model demonstrates the effect of 
aeroelastic deformation only. Each plot of the optimized models, thus, represents the effect of re-twisting of the wing 
and the aeroelastic deformation due to the re- twisting. 

X. Conclusion 

This study presents the development of a static aeroelastic model of a flexible sub-scale wind tunnel model based 
upon scaling of a full-scale transport aircraft wing of the NASA GTM. A static aeroelastic framework is developed by 
coupling a structural model of the flexible wing and the aerodynamic model of the aircraft. The structural model is 
constructed using finite-element modeling of an equivalent one-dimensional simple beam model of the wing. Aero- 
dynamic modeling is conducted using a vortex-lattice solution. The resulting static aeroelastic model is a coupled 
finite-element vortex-lattice model capable of converging aeroelastic solutions for the flexible wing model and devel- 
oping flexible aircraft lift curves and drag polars. 

The static aeroelastic model is implemented on a full-scale wing model of the ESAC or GTM. In order to develop 
the sub-scale model, a scaling procedure is developed first by geometrically scaling the full-scale model to sub-scale, 
then by conducting aeroelastic scaling. Aeroelastic scaling is conducted to scale the torsional stiffness of the sub-scale 
model such that the relative twist to span ratio matches that of the full-scale model, and the bending stiffness is scaled 
such that the sub-scale model has 10% wing tip deflection. Additional aeroelastic scaling is conducted to take into 
account the lower dynamic pressure and mach number of the wind tunnel test relative to the full-scale model’s cruise 
flight condition. Heuristic scaling factors are added to the scaling to account for coupled bending-torsion effects due 
to aerodynamic stiffening/softening. A static divergence analysis is conducted on the final, fully scaled flexible wind 
tunnel model to evaluate the risk of static instability of the model in wind tunnel testing. 

A final design of the wind tunnel model is developed by re-twisting the model through optimization targeted at 
minimizing induced drag or maximizing L/D at the wind tunnel design test condition. The gradient-based optimization 
approach is based on utilizing one-dimensional line searches in conjugate search directions. Two optimizations are 
conducted; one where additional pre-twist is applied to the wind tunnel model by linearly interpolating between pre- 
twist values specified at discrete points along the wing corresponding to the wing trailing edge extensions break and 
tip, and one where additional pre-twist is applied based on a shape function inspired by Chebyshev polynomials. The 
resulting optimized pre-twist of the wing is able to reduce the drag coefficient at the design test condition by 6.444% 
when optimized by specifying twist at discrete points, and 6.506% when optimized by specifying twist through a 
shape function. 

The final result of the study is a flexible sub-scale wind tunnel model configuration with static aeroelastic similarity 
to the full-scale ESAC wing but with increased wing tip deflection and tailored for the design test condition. An 
aeroelastic model accompanies the developed configuration that can be used in future validation against wind tunnel 
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testing. A clean wing model is analyzed, but investigation of the VCCTEF control surface can further use the model 
developed in future studies. 
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